#include <sstream>
#include <iterator>
#include <vector>
#include <string>
#include <iostream>
#include <fstream>
#include <map>
#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
#include <math.h>
#include <time.h>
#include <string.h>
//#include <my_global.h>
#include "CLineFields.h"

bool removeB = true;
bool removeHeader = false;

// return true if parameters are right.
bool parseParameters (int argc, char *argv[])
{
	for (int i = 0; i < argc; i++) {
		if (strcmp(argv[i], "-i")==0) {
			removeB = false;
		} else if (strcmp(argv[i], "-h")==0) {
			removeHeader = true;
		}
	}
	return true;
}


int main(int argc, char *argv[])
{

	if (argc < 2) {
		cerr << "convert the sam file created from the stripped reference file to sam file created from unstripped reference file " << endl;
		cerr << "usage: " << argv[0] <<" result.sam  [options] > result.new.sam" << endl;
		cerr << "Options:" << endl;
		cerr << "\t-i:\tkeep 'b' in chr name" << endl;
		cerr << "\t-h:\tremove file header lines" << endl;
		exit(-1);
	}
	CLineFields fileSam;
	if (fileSam.openFile(argv[1])==false) {
		cerr << "Can not open sam file: " << argv[1] << endl;
		exit (-1);
	}
	parseParameters(argc, argv);
	// process sam file
	fileSam.readline(removeHeader);
	vector<string> locations;
	int readLength = 0;
	while (fileSam.endofFile()==false) {
		if (fileSam.field[0][0]=='@' || fileSam.field[0][0]=='#') {// it is the header
			if (removeB==false || fileSam.field[1].find("b.")==string::npos) {
				cout << fileSam.line << endl;
			}			
		} else {
			// convert  (chr3b.300, 2)  to (chr3, 301)
			GenomePosition gp = GenomeLib::toRef(fileSam.field[2], fileSam.field[3], removeB);
			if (readLength==0)
				readLength = fileSam.field[9].length();
			int mm = 0;
			if (fileSam.field.size()>=14)
				mm = atoi(fileSam.field[13].substr(5).c_str());
			if (gp.allele != 0) { // on non-ref alleles, maybe only partially cover the allele.
				int more_mis = abs(atoi(fileSam.field[3].c_str()))-readLength;
				if (more_mis>0)
					mm += more_mis;
			}
			locations.clear();
			if (fileSam.field[4]=="0" && fileSam.field[fileSam.field.size()-1].substr(0,5)=="XA:Z:") { 
				// mapping quality==0 and we have chance to fix it                
                //locations.clear();// clear the old locations first
                //if (fileSam.field[fileSam.field.size()-1].substr(0,5)=="XA:Z:")
    			fileSam.divide(fileSam.field[fileSam.field.size()-1].substr(5),";",locations);

				for (int i=0; i<(int)locations.size(); i++) {
					if(locations[i].length()==0) {
						locations.erase(locations.begin()+i);
						i=-1;//invalidate the for loop.
						continue;
					}
					vector<string> f;
					fileSam.divide(locations[i],",",f);
					GenomePosition gp2 = GenomeLib::toRef(f[0], f[1], removeB);
					int new_mis = atoi(f[3].c_str());
					if (gp2.allele != 0) { // on non-ref alleles, maybe only partially cover the allele.
						int more_mis = abs(atoi(f[1].c_str()))-readLength;
						if (more_mis>0)
							new_mis += more_mis;
					}
					if ((gp.chrName==gp2.chrName && gp.pos==gp2.pos) || new_mis != mm) {
						// eg. the read is mapped to both (chr1, 201) and (chr1b.200, 2).
						// or the mismatch != the first hit
						if (new_mis < mm) {
							// this hit is better than the first hit, use this to replace the first hit.
							gp.chrName=gp2.chrName;
							gp.pos=gp2.pos;
							gp.allele=gp2.allele;
							mm=new_mis;
						}
						locations.erase(locations.begin()+i);
						i=-1;//invalidate the for loop.
					}						
				}
				if (locations.empty())
					fileSam.field[4]="23"; // make up the mapping quality
			}
			cout << fileSam.field[0] << "\t" << fileSam.field[1] << "\t" << gp.chrName << "\t" << gp.pos;
			for (int i=4, size=fileSam.field.size(); i<size; i++) {
				cout << "\t" << fileSam.field[i];
			}
			cout << endl;
		}
		fileSam.readline(removeHeader);
	}
	fileSam.fp->close();
}


